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This PhD tutorial is concerned with a description of the two-dimensional generalized Lorenz-Mie 
theory (2D-GLMT), a well-established numerical method used to compute the interaction of light 
with arrays of cylindrical scatterers. This theory is based on the method of separation of variables 
and the application of an addition theorem for cylindrical functions. The purpose of this tutorial is 
to assemble the practical tools necessary to implement the 2D-GLMT method for the computation 
of scattering by passive scatterers or of resonances in optically active media. The first part contains 
a derivation of the vector and scalar Helmholtz equations for 2D geometries, starting from Maxwell’s 
equations. Optically active media are included in 2D-GLMT using a recent stationary formulation 
of the Maxwell-Bloch equations called steady-state ab initio laser theory (SALT), which introduces 
new classes of solutions useful for resonance computations. Following these preliminaries, a detailed 
description of 2D-GLMT is presented. The emphasis is placed on the derivation of beam-shape 
coefficients for scattering computations, as well as the computation of resonant modes using a com¬ 
bination of 2D-GLMT and SALT. The final section contains several numerical examples illustrating 
the full potential of 2D-GLMT for scattering and resonance computations. These examples, drawn 
from the literature, include the design of integrated polarization filters and the computation of 
optical modes of photonic crystal cavities and random lasers. 
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I. INTRODUCTION 

Controlling the flow of light on the micro- and nano¬ 
scale level is one of the main technological breakthroughs 
enabled by the use of complex photonic media [1] . These 
complex media can be viewed as composite optical ma¬ 
terials whose characteristic length scale is of the order of 
the operating wavelength. Some prime examples of com¬ 
plex photonic media are photonic crystals mm , metama¬ 
terials [1] as well as random lasers mm- Besides the fact 
that the feasibility of all these examples has been demon¬ 
strated experimentally, a large body of work concerning 
the numerical modelling of complex photonic media has 
accumulated over the years. For several reasons, the de¬ 
velopment of these numerical tools is a challenging re¬ 
search topic in itself. Firstly, most of the tools typically 
taught in graduate courses in optics, such as the paraxial 
approximation and ABCD matrix approaches, are of¬ 
ten incompatible with wavelength-scale structures. Ad¬ 
ditionally, novel photonic media often involve a plethora 
of physical phenomena which can only be addressed nu¬ 
merically. A non-exhaustive list includes bandgaps [iia, 



2 


resonances [7] , slow light [S] , negative refraction [5] , topo¬ 
logical insulation [10] and Anderson localization of light 
m- Accordingly, numerical modelling tools aimed at 
complex photonic media must often compose with the full 
wave nature of the electromagnetic (EM) held (described 
by Maxwell’s equations), potentially active or non-linear 
media, as well as non-trivial boundary conditions of sub¬ 
wavelength geometries. 

Numerical schemes for the solution of Maxwell’s equa¬ 
tions that satisfy the above requirements can be classihed 
in two broad categories. The hrst category of methods 
is based on spatial and temporal discretization of the 
underlying differential equations m ESI- Well known 
examples include hnite element methods (FEM) as well 
as hnite-difference-time-domain approaches (FDTD). Al¬ 
though these fully numerical schemes are very general, 
meshing is most often associated with heavy memory re¬ 
quirements, especially in the case where the spatial ex¬ 
tent of the complex geometry considered is several times 
larger than the operating wavelength. In contrast, the 
second category of methods is based on the expansion 
of the EM field over a set of basis functions suited to 
specific geometries. In other words, these methods ex¬ 
ploit an underlying geometrical symmetry of the prob¬ 
lem. One obvious advantage of using a natural basis of 
the problem is that it alleviates the numerical difficulties 
associated with the first category, since it does not imply 
a spatial discretization of the problem space. The main 
topic of this tutorial is a method of the second category, 
which we refer to as two-dimensional generalized Lorenz- 
Mie theory (2D-GLMT) [T4|. This method is best suited 
for the modelling of effectively 2D complex photonic me¬ 
dia based on arrangements of disconnected and parallel 
cylindrical scatterers. These structures include, for in¬ 
stance, photonic devices based on coupled scatterers in 
the form of cylindrical holes in a dielectric slab waveg¬ 
uide, as well as devices based on coupled pillars m or 
coupled circular cavities [161 EZj- These effectively 2D 
devices are of great practical interest since they can eas¬ 
ily be integrated on photonic circuits. This allows the 
fabrication of innovative optical elements onto existing 
technological platforms, for instance silicon-on-insulator 
wafers [IIIHI or plastic substrates m- 

The purpose of this tutorial is to give a pedagogical, 
hands-on introduction to the theory behind 2D-GLMT, 
the main numerical aspects of the method as well as its 
wide range of applications. It is intended for graduate 
students or researchers wishing to construct their own 
implementation of the method and apply it to their re¬ 
search projects. The applications targeted in this tuto¬ 
rial are varied, but the emphasis is placed on two classes 
of optical geometries. The first is photonic-crystal in¬ 
spired devices, designed with the conventional wisdom 
of 2D photonic crystals (PhGs) in mind, but which are 
not necessarily infinite nor periodic. This research field 
has been active since the seminal 1987 publication by 
Yablonovitch, in which the author observed that multi¬ 
ple scattering and interference processes in a PhG could 


slow down the group velocity of light to the point of com¬ 
pletely inhibiting propagation |20| . Frequency intervals 
in which propagation is inhibited are called bandgaps, 
following the analogy with electronic bandgaps in semi¬ 
conductor crystals. Nowadays, PhG-inspired optical de¬ 
vices can readily be engineered using micro-fabrication 
methods such as UV m, holographic [H] or electron 
beam mm lithography. This relative ease of fabri¬ 
cation has enabled the realization of various integrated 
optical elements exploiting bandgaps in PhCs, includ¬ 
ing Mach-Zehnder interferometers [H[ |53|, waveguides 
pn E5] , mirrors [SSj , lenses [HllTj, polarization beam¬ 
splitters [5S], beam shapers [2S] and all-optical switches 
[29]. 

The second class of applications is called complex las¬ 
ing media, and include PhC lasers ED], photonic atoms 
and molecules [311132] as well as random lasers [61131- 
[55] . Although PhG lasers could fall in the first class 
of applications, we make the distinction since there are 
subtle differences between using 2D-GLMT for modelling 
passive devices and its use to compute lasing modes. 
Clearly these two classes do not exhaust the list of possi¬ 
ble applications of 2D-GLMT. Indeed, the method could 
be also extended to any optical geometry constituted 
of coupled homogeneous cylinders, for instance micro- 
structured hollow-core fibres [551155] . 

The tutorial is organized as follows. In section [ll] we 
review the basic concepts of electromagnetism relevant to 
2D-GLMT. We begin by deriving the vector and scalar 
Helmholtz equations for coupled cylindrical geometries, 
starting from Maxwell’s equations. After a derivation of 
the Helmholtz equations for passive optical media, we 
present a description of complex lasing media based on 
a stationary formulation of the Maxwell-Bloch equations 
called steady-state ab initio laser theory (SALT). This 
recently formulated approach, as well as the new kind of 
optical eigenmodes it introduces, are described in section 

imi 


Following these derivations of the basic electromag¬ 
netic equations for passive and active cylinders, we move 
on to a detailed description 2D-GLMT in section [TlT] This 
includes the derivation of expansion coefficients for two 
incident excitations, namely a plane wave (section ] 
and a collimated quasi-Gaussian beam (section 


HIEl) 


HE2I. 


The types of computations that can be achieved using 
2D-GLMT fall in two broad categories, which we refer 
to as “scattering” and “eigenmodes”, respectively. Since 
the objective is to describe as precisely as possible how 
to construct a numerical implementation of the method, 
examples of both types of computations are discussed in 
section m Specifically, PhC-inspired devices dedicated 
to beam shaping and polarization filtering are presented 
in section jlV A] as a representative example of scattering 
computations. Photonic c rystal cavities (section IV B) 
and random lasers (section IV C) serve thereafter as ex¬ 
amples of eigenmode computations. 
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II. ELECTROMAGNETIC THEORY 


in (4a) and (4b) yields, after a few vector identities, 


Our starting point is Maxwell’s equations which, for 
a non-magnetic medium (unit magnetic permeability, or 
fi = 1), take the form 


V • D = p, 

V • H = 0, 

1 aH 

c dt ’ 

19D J 

c dt c' 


V X E = - 
V X H = 


(la) 

(lb) 

(l c) 

(l d) 


Without loss of generality, the system of Heaviside- 
Lorentz electromagnetic units is used throughout this 
work except where indicated. Details of the conversion 
between Heaviside-Lorentz and other systems, i.e. SI or 
Gaussian units, can be found in |39j . 

We are mainly concerned with obtaining rigorous so¬ 
lutions of Maxwell’s equations in regions of space that 
contain no free charges and no free currents, that is [40] 


f fP‘'W 

V(V.E)-V^E+-—=0, (6a) 

V(V.H)-V^H+A^-^|[ExVx.]=0, 

(6b) 


where e = 1 -I- Xe denotes the relative dielectric permit¬ 
tivity. When e(r) is a continuous function of space, these 
equations can not be further simplified. However, in 
the case of arrays of individually homogeneous cylinders, 
e(r) is a piecewise constant function. In other words, 
we will be concerned with regions of constant relative 
permittivity separated by well-defined interfaces, where 
the value of e(r) is discontinuous. It is then possible to 
solve the equations under the assumption that Ve(r) = 0 
and Vxe(r) = 0, and connect the solutions at the cylin¬ 
ders’ interfaces using appropriate electromagnetic bound¬ 
ary conditions. Under these conditions, (la) reads 


V • D = V • [eE] = eV • E -f E • Ve = 0 


(7) 


p = 0, (2a) 

J = 0. (2b) 

However, since we are also interested in solutions of 
Maxwell’s equations in active media, we suppose the fol¬ 
lowing relation between the displacement field D and the 
electric held E [311 HO] 


D = E -|- Pn 


(3) 


where Pt is the total polarization held. It contains both 
a contribution from the linear response of the medium 
(basically the optical density, or refractive index) and the 
non-linear response. Taking the curl of ( [I^ and (Id), to¬ 
gether with ^ , one obtains the following wave equations 


V X V X E 


1 a^E 1 a^p 

+ —= 0 , 


at2 


dt^ 


Vx VxH+4^--^[VxPt] =0. 

ot^ c at'- ■' 


(4a) 

(4b) 


A. Dielectric media 


The two fundamental wave equations for E and H can 
be further reduced if the response of the optical media 
is assumed to be linear. In this case, the polarization 
held Pt is a linear function of the electric held, more 
specihcally 


Pt = XeE, (5) 

where Xe is the electric susceptibility of the medium, i.e. 
Xe = 0 in free space. The substitution of the Ansatz ([^ 


and since Ve(r) = 0, one has immediately that V • E = 0. 
This simplihes (6a I and (6b) to the vector Helmholtz 
equations 


r, e(r) 

a^E 

(8a) 

V^E 

dt^ 

o efr) 


(8b) 

V^H 

dt 


These equations are more conveniently solved in the fre¬ 
quency domain, using the following substitution for tem¬ 
poral derivatives 


dt'^ 


i - UJ . 


(9) 


This is of course the assumption of a harmonic time- 
dependence of the electromagnetic helds, i.e. E = 
E(r)e-*‘^‘ and H = H(r)e-*‘^‘. This last step yields 

V^E-f e(r)fc2E = 0, (10a) 

V2H + e(r)fc2H = 0, (10b) 


where the wavenumber k = uj/c. 

The vector Helmholtz equation is a general description 
of wave propagation in fully three-dimensional, linear me¬ 
dia of piecewise constant refractive index (Ve(r) = 0). 
This covers a wide range of optical geometries: opti¬ 
cal hbres and waveguides m, dielectric cavities such 
as spheres and tori [32| and three-dimensional PhCs 
[SlElllis] to name a few. 

Since the scope of this tutorial is cylindrical geome¬ 
tries, an additional symmetry can be used to further 
reduce the vector Helmholtz equation to a scalar form. 
Effectively two dimensional geometries suggest two priv¬ 
ileged polarization directions. Formally, it is possible to 
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uncouple the electromagnetic field in two orthogonal po¬ 
larization components if variations of the relative permit¬ 
tivity e(r) are restricted to a plane, say the (x^y) plane, 
meaning 


e(r) = e(x,y). 


( 11 ) 


This 2D hypothesis implies that the electromagnetic field 
can be uncoupled in a transverse-electric (TE) compo¬ 
nent (E field component parallel to the a;, y plane, H field 
component normal to the x^y plane) and a transverse 
magnetic (TM) component (H held component parallel 
to the X, y plane, E held component normal to the x, y 
plane). Under this two-dimensional restriction, (10aI and 
(10b I become a single scalar Helmholtz equation, that is 


V'^y}{x, y) -h e{x, y)k‘^(p{x, y) = 0, 


( 12 ) 


where (p stands for either the or held component. 
For the remainder, we shall also suppose that optical ge¬ 
ometries have inhnite extent in the direction perpendic¬ 
ular to the (x, y) plane and no propagating component 
in the z direction. This allows one to ignore the z de¬ 
pendence of the EM held. We shall also assume that op¬ 
tically thin geometries, such as slab waveguides, can be 
taken into account via an effective value of the permit¬ 
tivity e(x, y). As with (10a) and (10b), the scalar form of 


the equation means the eigenfunctions of the Helmholtz 
equation are of similar form for both polarizations, with 
however differing boundary conditions (see for instance 
section |III C). Let us recall once again that equations 


(10aI, (10b) and (12) apply only for constant or piece¬ 


wise constant refractive index. 


1. Scattering vs resonances: Quasi-bound states 


The linear Helmholtz equation (12) may be used for 
two kinds of computations, namely wave scattering and 
resonances. As will be discussed in section [HI E[ scatter¬ 
ing computations using 2D-GLMT are straightforward 
and consist in computing the response of a given opti¬ 
cal system to an incident excitation. The incident and 
scattered wave are expanded on a suitable basis of the 
Helmholtz equation (in the case of 2D-GLMT, cylindri¬ 
cal harmonics) and the expansion coefficients of the scat¬ 
tered wave are deduced from the expansion coefficients of 
the incident excitation. Thus, wave scattering problems 
imply solving a inhomogeneous linear system of equa¬ 
tions. 

On the other hand, resonance computations are more 
subtle since they are analogous to solving a homoge¬ 
neous system, in other words finding the non-linear eigen¬ 
values of a system of equations. The associated solu¬ 
tions of the Helmholtz equation are consequently termed 
eigenmodes. Gomputing the eigenmodes usually implies 
searching for complex k values {k = k' -\- ik") satisfying 
a resonance condition. Complex k values with k” < 0 
result in exponential decay of the total energy of the 


optical system. According to the harmonic hypothesis 
p{v,t) = (/3(r)e“*“*, the total electromagnetic energy of 
the system is proportional to 

\p{v,t)\^ = W{v)\^e-^/\ (13) 

where |</9(r)P is the field profile of the eigenmode and r 
the characteristic decay time (or half-life) of the system 
energy. It is given by 


2ck" 


(14) 


This delay can be readily interpreted as the characteris¬ 
tic residency time of a photon trapped in an optical res¬ 
onator. Since these modes are inherently leaky, they are 
usually called meta-stable or quasi-bound (QB) states, by 
opposition to bound states characteristic of closed quan¬ 
tum systems. This analogy stems from the similarity be¬ 
tween the Schrodinger equation and the Helmholtz equa¬ 
tion [iij . 

Although useful for modelling active cavities, QB 
states come with an important difficulty. Whereas the 
solutions of (12) are assumed to be stationary inside the 


cavity, resonant wavefunctions of open systems cannot 
by definition be stationary since their energy must decay 
in time |45j . Since e(r) is real everywhere, the only way 
to obtain a decaying energy is through the imaginary 
part of the eigenfrequency k. However, this imaginary 
part has the net effect of introducing an artificial gain in 
the whole space domain. Consequently, QB states grow 
exponentially towards inhnity, which is not physically re¬ 
alistic [15]. Strictly speaking, these states are not even 
regular functions because of this blow-up behaviour [1^ . 
In section HB[ we introduce a recently formulated theory 
that overcomes this shortcoming by taking into account 
the active medium ab initio, instead of the a posteriori 
introduction of gain associated with QB states [15] . 


B. Optically active media: Steady-state ab initio 
laser theory (SALT) 


In the previous section, we have derived the scalar 
Helmholtz equation (12) for 2D geometries. The main 


hypothesis behind the derivation is the linear response of 
the optical media. In the presence of an active medium, 
for instance a pumped laser cavity, the derivation is some¬ 
what different. As discussed in section [HA 1[ it is still 
possible however to use the linear Helmholtz equation 
to model 2D active media using QB states. This de¬ 
scription comes with two main shortcomings: the ex¬ 
ponential growth of QB states outside the resonator as 
mentioned previously and the conspicuous absence of the 
gain medium parameters in the linear Helmholtz equa¬ 
tion ([T^. 


To circumvent this difficulty, we will make use of a re¬ 
cent formulation called steady-state ab initio laser theory 
(SALT). The term ab initio refers to the fact that SALT 
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only requires the distribution of e(r) of the passive cav¬ 
ity, or resonator, and a number of parameters describing 
the gain medium. The theory is stationary, meaning it 
works in the frequency domain, and is intended to bridge 
the gap between the over-simplified QB states approach 
described in section III A II and time-domain simulations 
using dynamical theories, for instance the Maxwell-Bloch 
or Schrodinger-Bloch theories [3Sl US]- Partly equiva¬ 
lent formulations devised for the extraction of the lasing 
thresholds and frequencies of optical cavities are also pre¬ 
sented in Refs. 

SALT was initially developed by H. Tiireci, A. Stone 
and B. Collier in 2006, and originally called ab initio self- 
consistent laser theory, or AISC |46j . It was this seminal 
paper that first introduced the constant-flux (CF) states 
central to the ab initio approach. The theory was fur¬ 
ther developed in [34l [50U52] and renamed SALT in the 
process. Further generalization of SALT to complex gain 
media is to this day an active research topic . In sec¬ 
tion IIB 1 we highlight the defining features of the theory 
and present the governing equations for the eigenmodes 
of a 2D laser cavity near threshold. Since the CF states of 
an array of active cylinders can be readily computed us¬ 


ing 2D-GLMT, we discuss them further in section IIB 2 


fields in a positive frequency part and a negative fre¬ 
quency part, i.e. E = -\- E~ and P = P'^ -\- P~. This 

allows to rewrite ( [T7| ) only for the positive frequency field 
component, keeping in mind that the negative frequency 
component satisfies a similar equation 


V'^E+ 


e d^E+ 

(? dP 


1 d^P+ 

dP 


= 0 . 


(18) 


The SALT theory is based on the Maxwell-Bloch equa¬ 
tions for a two-level atomic system |46] : 

gp+ (j2 

— = -{loj, + 7j.)P+ + ^E+V, (19) 

^ = 7||(2?o(r) -V)- ^{E+iP^r + P+iE+n 

( 20 ) 


These equations relate the electric and polarization fields 
with the spatially varying population inversion 'D{r) and 
the pump profile T’o(r). The gain central frequency uja 
and linewidth (or polarization relaxation rate) 7 _l appear 
in (19). Other constants are the population relaxation 
rate yy and the dipole matrix element g. For convenience, 
we renormalize the population inversion 


1. Basic equations and threshold lasing modes 

In this section, we explain how to obtain the basic 
equations of SALT using the electromagnetic theory pre¬ 
sented in section |ll] and the Maxwell-Bloch equations for 
a two-level atomic system. Since the goal is not to present 
a detailed discussion of SALT but rather to describe the 
eigenstates that are compatible with 2D-GLMT, we shall 
only present an abridged derivation of the SALT equa¬ 
tions and briefly discuss the physical properties of the 
eigenstates. For convenience, only the derivation for a 
TM polarized wave and a 2D optical geometry is pre¬ 
sented. We refer the interested reader to [35] or |3I| for 
a complete theoretical description of SALT. 

The first step in the derivation of the SALT equations 
is to introduce a non-linear polarization term of the fol¬ 
lowing form 


Pt = XeE -I- P 


NL5 


(15) 


in (4a), which yields 
V X V X E - 


e 1 


NL 


dP 


dP 


= 0 . 


(16) 


Under the assumption of a piecewise constant refractive 
index, and for TM polarization with E = E^, (16) sim¬ 
plifies to 


e d'^E 1 d^P 

P dP P dP 


(17) 


where P = {Pnifjz- Under the rotating wave approxi¬ 
mation [52] , one can expand the electric and polarization 



( 21 ) 


so that D is dimensionless. This results in the following 
Maxwell-Bloch equations 


dp+ 

—— = -{iuJa + 1A.)P^ - ( 22 ) 

ot 

= {Do{r)-D) + iK\E+{P+r ^P+{E+r], 

7|| dt 

(23) 


where the coupling constant n is defined as 

7||7_lS2' 


(24) 


This coupling coefficient governs the strength of the 
spatial-hole burning effect, that is interaction between 
lasing modes above threshold. 

The main hypothesis of SALT resides in the so- 
called stationary inversion approximation, which allows 
to recast the Maxwell-Bloch equations in a numerically 
tractable form. In short, this approximation implies that 
the population inversion of the two-level system remains 
constant in time (^ = 0). This approximation is valid 
in the single-mode regime, as well as in the multi-mode 
regime under the condition 7 j_ 3> yy, if the character¬ 
istic time scale of the inversion dynamics is larger than 
that of the polarization dynamics. Furthermore, in the 
single-mode regime or when the electromagnetic fields 
are small with respect to the coupling coefficient k, the 
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hole-burning effects can be neglected. In this case, a mul- 
tiperiodic Ansatz yields a set of threshold lasing modes, 
or TLMs: 




—icku.t 


— ickn,t 


(25a) 

(25b) 


This set of modes is governed by the following equation 




e(r) 


laPUFir) 
kf,-ka+ i-fa. 


(^^(r) = 0 , (26) 


where ka = Wa/c is the gain centre frequency (we use fre¬ 
quency and wavenumber interchangeably), 7 ^ = 7 _l/c is 
the gain width, Dq is associated to the “pump strength” 
and F{r) is the spatial pump profile. In the absence 
of pumping (F(r) = 0 ), for instance outside the laser 
resonator, the modes satisfy a linear Helmholtz equa¬ 


tion (12 1 with real frequency k. This is a more realistic 


description of the lasing modes than the QB states ap¬ 
proach since the gain medium is introduced as an effective 
complex-valued permittivity inside the resonator instead 
of an artificial imaginary part of the wavenumber |46] . 

The determination of TLMs can be achieved using 
various numerical schemes including the Lorenz-Mie ap¬ 
proach described in section |III| or the finite element 
method described in [54j . For a given combination of 
e(r), F(r), ka and 7 a, the problem consists in finding 
the value of the lasing frequencies thresholds Dq and 
the field distributions However, for a given choice of 
exterior real frequency k, the value of Dq is in general 
complex. Consequently, the TLM must satisfy an addi¬ 
tional reality condition on the threshold value, since a 
complex-valued Ug does not correspond to a physically 
realistic lasing mode. This means that one must map 
the evolution of Dq as a function of k, and when the 
former crosses the real axis at k = k^, obtain a pair of 
real numbers (fc^, Dq) defining the TLM lasing frequency 
and lasing threshold, respectively. The first lasing mode 
is therefore the TLM with the smallest threshold Dq [52] . 


2. Constant-flux states 

One of the defining features of SALT is the introduc¬ 
tion of new eigenstates, such as the TLMs described by 
(26). The theory also introduces another kind of eigen¬ 
state called a constant-flux (CF) state jUj. CF states 
are useful because they are straightforward to compute 
and can be used as a basis to expand TLMs [S2]. They 
are parametrized by real wavenumbers outside the laser 
resonator, and are thus physically meaningful when com¬ 
pared to QB states. The basis of CF states satisfy the 
following modified Helmholtz equation 


[V2 + e{r)Kl{k)]g,^ = 0 


r G C, 


[V2 + eir)k^]ip^ = 0, r^C, 


(27a) 

(27b) 


where C is the cavity region, defined as the union of all 
optically active regions. The eigenvalues are complex 
and depend on the exterior frequency k, which is always 
real. This formulation ensures that the total electromag¬ 
netic flux outside the cavity is conserved [^ . 

In the special case of a uniform pumping inside the 
cavity region, i.e. F(r) = 1 for r G C and F{r) = 0 
for V C, and a uniform refractive index distribution 
e(r) = Cc for r G C, the TLMs can be trivially exp anded 
in the basis of CF states. By comparing (26) and (27a), 
one obtains the relation 




k-kaF i^a 


- 1 


(28) 


In other terms, if the active medium is uniformly 
pumped, one can simply compute the CF states of the 
geometry, given the value of Cc, and extract the complex 


values Dq of the associated TLMs using (28). Moreover, 


in the simplest possible approximation, CF states may 
be regarded as the actual lasing modes, for instance in 
the case of a gain transition frequency close to the CF 
state eigenfrequency [15] . 


C. Circular cavity: QB vs CF states 

This section completes the description of QB and CF 
states by giving their explicit expressions for the sim¬ 
plest 2D geometry for which the solution is analytic: a 
single homogeneous dielectric cylinder, or circular cav¬ 
ity. Since this tutorial is mostly concerned with coupled 
cylinder geometries, it is also useful to present the closed- 
form expressions for the eigenmodes of a single, uncou¬ 
pled cylinder, which is a limiting case of the Lorenz-Mie 
theory presented in section jlH] 

Consider a dielectric circle of refractive index Uc = 
^Jef and radius r placed in a medium of refractive index 
rig. The cavity is centred at the origin of a “global” 
cylindrical coordinate system. Under these conditions, 


the QB states of the circular cavity, solutions of Eq. (12), 
are 111 ] 




AiJi{nckp), 


BiH, 


(+) 


(nokp). 


p<r, 
p> r, 


(29) 


where I is the angular quantum number of the solution, 
Ji is a Bessel function of the first kind and = Ji-\-iYi 
is a Hankel function of the first kind. Analogously, the 
CF states of the circular cavity, solutions of Eq. (27), are 


f Q\ I AiJi{naK{k)p), p<r, 

[BiH) ’{nokp), p>r. 


(30) 


One notes that the CF eigenvalue K{k) explicitly appears 
in the solution for p < r. 
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There exists an infinite set of QB and CF states for 
any cavity. In the case of the circular cavity, both types 
of states can be characterized by a pair of quantum num¬ 
bers {l,j), corresponding to the number of angular and 
radial lobes of the eigenmode, respectively. Moreover, the 
expressions given by Eqs. (29) and (30) satisfy an outgo¬ 


ing boundary condition outside the resonator, sometimes 
called Sommerfeld radiation condition. This condition 
ensures a net flow of the electromagnetic energy of the 
eigenstate towards infinity, and is mathematically given 
by 


„ikp 

lim ^{p,e) = 
p^oo \/kp 


(31) 


In order to determine the eigenvalues of QB states, one 
must find non-trivial solutions to a homogeneous linear 
system relating the unknown coefficients Ai and Bi. This 
linear system, obtained via the application of electromag¬ 
netic boundary conditions (see section III C) at p = r is 
given by 


Jiiuckr) -Hl~^\nokr) Ai^ 

c„UcJ/(ncfcr) -<roUoij/^^'(nofcr) 



= 0, (32) 


where = 1 for a TM (TE) polarized wave. 

These factors account for the polarization of the mode. In 
essence, one must find the complex eigenvalues k = ^qb 
such that the determinant of the coefficient matrix in Eq. 
(32) is zero. Once this is done, the Ai and Bi coefficients 


are readily obtained from the linear system. A similar 
system of equations can be obtained for the CE states of 
the circular cavity 


/ JiiricKr) -H['^\nokr) \ /AA 

WnncKJ'i{ncKr) -<;onokHj;~'~'’'(nokr)j \Bi) 


= 0 . 


(33) 

In this case, one must hnd eigenvalues K^{k) for a fixed 
exterior frequency k. 

Eor comparison, a QB state of the circular cavity and 
the associated CE state are plotted in Fig. Both 
modes are characterized by the pair of quantum numbers 
{l,j) = (10,3). Since we chose the exterior frequency of 
the CF state close to the real part k'q^ of the QB eigen¬ 
value, the field profile inside the cavity (p < r) is similar 
for both eigenstates. However, one clearly sees that the 
CF state remains bounded outside the cavity, whereas 
the QB state begins to grow exponentially roughly be¬ 
yond p = 1.5r. This behaviour is consistent with the 
Sommerfeld radiation condition. If one injects a complex 
k value in Eq. (31) with k” < 0, the solution clearly 


blows up at infinity. In contrast, CF states are always 
characterized by real wavenumbers outside the cavity re¬ 
gion, which ensures that the solution remains bounded. 
This example of the circular cavity serves two purposes. 
Firstly, it shows that QB and CF states can be uniquely 
mapped onto each other for the corresponding quantum 
numbers. Actually, this mapping is also possible in the 


case of more complex cavity geometries |46j , a prime ex¬ 
ample of which is random lasers, discussed for instance in 
[55H55] . Secondly, it highlights the unphysical behaviour 
of the usual QB states at infinity, a behaviour which can 
be corrected using the formulation of CF states. 

To conclude this section on SALT, the main features 
of the three kinds of eigenstates (QB states, CF states 
and TLMs) are summarised in Table |l] 


III. 2D GENERALIZED LORENZ-MIE THEORY 

We now describe in some details the theoretical frame¬ 
work which we refer to as 2D-GLMT. According to 
Gouesbet and Lock El, the expression “generalized 
Lorenz-Mie theories” (GLMTs) is used to generically 
denote a set of light-scattering theories describing the 
interaction between an incident electromagnetic beam 
and a discrete scattering particle (or array of parti¬ 
cles) provided that the scatterers possess enough sym¬ 
metry, allowing to solve the problem using the method 
of separation of variables. The approach has been re¬ 
ferred to by a number of names in the past, including 
“multipole method” or “fast multipole method” in Refs. 
[33l|36l|37l[56]. Historically, one should mention that nei¬ 
ther Ludvig Lorenz nor Gustav Mie were involved in the 
solution of Maxwell’s equations of coupled cylinders. The 
solution of Maxwell’s equations in cylindrical geometries 
is best associated to Lord Rayleigh m, and some au¬ 
thors prefer the term “Rayleigh scattering” for cylinders 
instead of “Lorenz-Mie scattering” or “Mie scattering” 
which usually deals with spheres [14]. We will avoid the 
term “Rayleigh scattering” since it most often refers to 
scattering by particles much smaller than the wavelength 
of the incident light, an approximation not used in the 
present treatment. 

This section presents a “reference implementation” of 
the 2D-GLMT and its application to electromagnetic 
modelling of two-dimensional arrays of cylinders. The 
derivations and notation used follow references musa- 
I60j . Similar derivations are also found in a variety of 
other published works including 1331111. 

This section is organized as follows. The first sub¬ 
sections (sections HI A to HIC) are concerned with 


the basic theoretical treatment of 2D-GLMT. In section 
|HI A[ we state Graf’s addition theorem for cylindrical 
functions since it is at the heart of 2D-GLMT. Section 
|HIB I presents the basic equations and the main hypoth¬ 
esis of 2D-GLMT, the expansion of the electromagnetic 
fields in a basis of cylindrical waves. To isolate the un¬ 
known expansion coefficients, electromagnetic boundary 
conditions must then be enforced as described in section 

Imcl 

The next two sub-sections deal with two classes of 
computations enabled by the use of 2D-GLMT. In sec¬ 
tion |llT^ we describe scattering computations, for which 
2D-GLMT is an ideal tool to obtain the field produced 
by a finite array of cylindrical scatterers. We show how 



















Figure 1. Comparison of the field profile of a QB and a CF state of a circular cavity of refractive index ric = 1.5 embedded 
in a medium of refractive index no = 1. (TM polarization, arbitrary normalization). The eigenmodes are characterized 
by the quantum numbers {l,j) = (10,3). The QB state eigenvalue is fege = 13.521 — 0.442z. The CF state eigenvalue is 
Ki_i = 13.558 — 0.440i (exterior frequency k = 13.52). Reproduced from [55]. 


Table I. Comparison of the three kinds of eigenstates used in this tutorial. Outside the cavity region, all eigenstates are governed 
by the linear Helmholtz equation (12l. The defining features listed in this table assume that the refractive index distribution 
e(r) is a real and piecewise constant function. 


Name Acronym Governing equation inside cavity region Defining features 


Quasi-bound state 

QB 

[V2-te(r)fc2]vJ^(r) = 0 

Grows exponentially at infinity; k 
complex 

Constant-flux state 

CF 

[V^ + e{r)Klik)]Mr) = 0 

Remains bounded at infinity; 
complex, k real 

Threshold lasing mode 

TLM 

|V2 + 

,(r)+ ■ 

ka 

f~ 0 Includes a description of the pump 
' and gain medium; real 


this formalism is compatible not only with the scatter¬ 
ing of plane waves by cylinders, but also of non-paraxial 
Gaussian beams via the complex-source beam approach. 
Eigenmode calculations follow in section fill F[ where las¬ 
ing modes of an array of cylinders are presented. More 
specifically, we will be interested not only in the compu¬ 
tation of the usual quasi-bound (QB) states, but also of 
the constant-flux (CF) states, central to the SALT the¬ 
ory. An equivalent method for computing modes of ran¬ 
dom lasers has been summarily described in [33] under 
the name “multipole method”. 


A. Graf’s addition theorem 

We state Graf’s addition theorem [SH Eq. 9.1.79] since 
it is central to the derivation of the main equations of 
2D-GLMT. This addition formula allows one to displace 
Bessel functions from one cylindrical system of coordi¬ 
nates into another. Let J- denote J,Y, H^~'> or 
any linear combination of these functions. The follow¬ 


ing identity holds 

OO 

J-,+„([/)J™(E)e*™“, (34) 

m— — (x> 

where < jt/j and 

VF^ = {72 + i/2 _ 2 UVcosa, 

U — V cos a = W cos y, (35) 

V sin a = W siny. 

The branches must be chosen such that W —>■ U and 
X —f 0 as E —>-0. The restriction < \U\ is 

unnecessary if = J because the Bessel function of the 
first kind and integer order is a bounded entire function. 
If U, V, W are real numbers, they can be interpreted as 
edges of a triangle as shown in Fig. [2}r. The theorem 
holds for complex-valued arguments also. 
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Figure 2. Graf’s addition theorem, (a) Triangle interpreta¬ 
tion of the theorem, (b) Application of the theorem to Eq. 


(421. 


B. Basic equations 


fields can be expanded in a basis of cylindrical functions, 
that is 


V5-®(r) = vJ,(r)-hv5s(r), (37) 

with 

OO 

X! (38a) 

1 — — 00 
N OO 

(,Ss(r) = ^ ^ 6„pij/,’^^(fco/9n)e*'®". (38b) 

n—1 l' — — oo 

The coefficients {a°;} are the beam-shape coefficients, 
used to parametrize the incident field in the frame of 
reference of the scatterer. 

Inside the scatterer, the field is similarly repre¬ 
sented as 


Consider an array of N cylindrical scatterers of radii r„ 
and relative permittivity e„. Let also r„ = (p„, 0„) be the 
cylindrical coordinate system local to the scatterer, 
whose center is located at R„ = {Xn,Yn). We suppose 
that every cylinder (or hole) is infinite along the axial 
z direction. The field outside the cylinders satisfies the 
following Helmholtz equation 

+ kl]^{x,y) = t), (36a) 

with fco = whereas the field inside the n**' scatterer 

satisfies the following Helmholtz equation 

-f k1]ip{pn, On) =0, Pn< Tn- (36b) 

For generality, we do not suppose any relation between 
kn and e„ at this point. This will allow the derivations 
to be compatible both with passive and active media, as 
described in section nm 

The central hypothesis of 2D-GLMT is that the field 
exterior to the scatterers can be written as the su¬ 
perposition of an arbitrary incident beam and the sum 
of the field scattered by each individual scatterer. The 


OO 

ifnir) = ^ Ji(fc„/9„)e*'®". (39) 


In order to apply electromagnetic boundary conditions 
at the interface of the scatterer, one must find an ex¬ 
pression for containing only cylindrical harmonics 

centred on the scatterer, that is 


OO 

T’n ^ ^ (/coPn) 


bnlHj;~'~\koPn) 


1 — — 00 


(40) 

This can be achieved via the application of Graf’s ad¬ 
dition theorem for cylindrical functions, performing a 
translation from the frame of reference of scatterer n' 
to the frame of reference of scatterer n. We apply Graf’s 
addition theorem with an appropriate change of variables 
(see Fig. for the definition of angles and edges) 


U — k^Rjin'tY — kQPn,^Y — kt^Pn’• (^^) 


The theorem states that 




(42) 


1 — — 00 


where i?„„' is the center-to-center distance between scatterers n and n' and is the angular position of scatterer 
n' in the frame of reference of scatterer n. Substituting (42) in (38) yields 


OO OO 

p^{v)= Y a°iJ/(fcoPn)e*'®" + ^ bniHY{koPn)e^'^- 

1 — — 00 1 — — 00 

OO OO 

+ E E E Hl+l{koRnn')Jl{koPn)e^'‘>-. 

1——00 n'^n l'——oo 


( 43 ) 
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The comparison of (40) with (43) results in the following relation between the {a„/} and {bni} coefficients 


n'^n 


(44) 


r 


At this point, we should stress that 2D-GLMT is not 
restricted to disconnected cylinders embedded in an in¬ 
finite medium. Cylinders embedded in other cylinders 
are also within the reach of the method, provided appro¬ 
priate expressions of the field are used. The surround¬ 
ing cylinder may represent the jacket of an optical fibre, 
for instance, and the leaky modes of this geometry are 
readily computed using a slightly modified version of 2D- 
GLMT, as detailed in Refs. 157] . For the remainder, 
however, we shall restrict the discussion to disconnected 
cylinders with the field in the surrounding medium given 
by m. 


C. 


Electromagnetic boundary conditions and 
characteristic equation 


Relation (44) is one of the central results of 2D-GLMT. 


It accounts for the mutual influence of all individual scat- 
terers. In order to obtain the {uni}, {bni} and {cni} for 
a given set of beam-shape coefficients {a°;}, one must 
apply electromagnetic boundary conditions to (39) and 
([40| at pn = Tn- For a TM polarized wave, the first con¬ 


dition is the continuity of across the cylinder interface 
at Pn = Tn- The second condition is the continuity of the 
component of H parallel to the interface compo¬ 
nent). From ([I^, one obtains 


= -;^[V X E], 




For a TE polarized wave, the first condition is rather the 
continuity of across the cylinder interface at Pn = fm 
and the second condition is the continuity of the com¬ 
ponent of E parallel to the interface component). 

From (Id), one obtains 


^ i i dHz 


(46) 


The four conditions of continuity (two for each polariza¬ 
tion) can be rewritten as 






dpn 


= <^0 


Qpn 


(47a) 

(47b) 


Again, the q factors account for polarization. Specifi¬ 
cally, we have q = 1 (l/e^) for a TM (TE) polarized 
wave. Thus, applying this boundary condition to (39) 


and (40), one obtains 


T bnlH^ (^O^n); (48a) 

^nEn^nJli^n'^n) ~ *7/(^O^n) T ^nAO^O-^/ (^O^n)- 

(48b) 


Eliminating Cni, we obtain the relation bni = On/Sni, with 

Snliko,kn) = (49) 

where 


Hl+^'(korn)-rniHl+\korn)' 

kn^l {kn'^n) 


and 


bn/ — CnO' 


Cij — 1 


koJliknTn) ’ 


(50) 


(51) 


for a TM (TE) polarized wave. 


Substituting bni = aniSni in (44) yields 


bnl-SnlJ2 e<^'-^^^-’^'Hl^likoRnn')bn'l' =Sniali. 

n'^n l'——oo 

(52) 

This relation between the {a°;} and {bni} coefficients can 
be rewritten in matrix form as 


Tb = ao. 


(53) 


with 


TS = Snn'Sw - (1 - Snn')e^^''-'^^'^-’Hl^likoRnn')Snl, 

(54) 

and 


^0 — { SnlO^nl 


h=Hniy 


(55) 


The transfer matrix T is typically constructed by trun¬ 
cating the series expansions to order /max- It is composed 
oi N X N blocks of dimension 2Z„iax +1, where /max is the 
truncation order. To ensure an adequate representation 
of the field, the truncation order is often chosen as a linear 
multiple of the ratio between the characteristic dimension 
of the problem, for instance the radius of the larger cylin¬ 
der comprised in the array, and the wavelength. Irrespec¬ 
tive of the value of Imax, the size of the transfer matrix 
scales like where N is the number of cylinders con¬ 
sidered. This suggests that 2D-GLMT is best suited to 
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Figure 3. Example scatterer array consisting of three identical 
coupled cylinders with permittivity Cc = 4, arranged on the 
vertices of an equilateral triangle (-R 12 = 2.5ri). The transfer 
matrix of this arrangement is shown in Fig. Reproduced 
from [55]. 


cylinder dimensions of the order of the operating wave¬ 
length, in which case it provides fast and accurate results 
without the need for a spatial discretization of the prob¬ 
lem. Otherwise, in the case of cylinders much larger than 
the operating wavelength, and/or a very large number of 
scatters, the memory and processor requirements of the 
method may become impractical. In that case, alterna¬ 
tive methods may be more appropriate than 2D-GLMT, 
for instance the finite element method |54l |63] . 

As a representative example, the transfer matrix for 
an array of three identical (cc = 4) coupled dielectric 
cylinders arranged on the vertices of an equilateral tri¬ 
angle (Fig. is shown in Fig. This matrix is clearly 
composed of 3 x 3 blocks of dimension 2Zniax + 1; where 
the prescription Imax = 3fcri sets the truncation order 

[SHIES]. 


D. The characteristic equation as a Fredholm form 
of the second-kind 


Although 2D-GLMT in the form presented in the pre¬ 
vious section has been widely used for two decades, it 
comes with an important shortcoming. Specifically, the 
underlying matrix equation (531 cannot be truncated at 


an arbitrary order to provide convergence. In this sec¬ 
tion, we isolate the source of this difficulty and show how 
it can be completely alleviated using a straightforward 


procedure. This has been acknowledged in recent publi¬ 
cations including [64], and generally solved in [61] . 

As can be demonstrated by using the asymptotic form 
of cylindrical functions for large positive order [HI], for 
large values of |F|, the elements of the off-diagonal blocks 
of grow exponentially. On the one hand, for large 
values of \l\, the elements of the off-diagonal blocks of 
decay exponentially. The net result is that, for 
small values of \l — l'\, the elements of the off-diagonal 
blocks remain bounded, as the exponential growth with 
respect to |Z'| is compensated by the exponential decay 
with respect to |Z|. On the other hand, for large values of 
|Z —/'I and large jl'l, the matrix elements quickly blow up. 
This exponential growth can be clearly seen in Fig. & 
for off-diagonal elements. The consequences are impor¬ 
tant, since it is theoretically not possible to minimize the 
error of 2D-GLMT by solving arbitrarily large matrices 
because it would involve arbitrarily large numbers. 

To avoid numerical instabilities, one can follow the 
rule of using “not too large” truncation orders. For a 
given cylinder array, the most common prescription for 
the truncation order Z^ax is [SHI El] 

^max 3A:rniaxi (56) 


where r^ax is the radius of the largest cylinder in the ar¬ 
ray. With this rule, it is usually possible to solve scatter¬ 
ing and resonances problems with an acceptable degree 
of accuracy. 

Despite the fact that the prescription of using “not 
too large” truncation orders is adequate in many in¬ 
stances (see Ref. [SH] EHl ISI] and references therein), a 
robust method guaranteeing the convergence of the ma¬ 
trix equations to an arbitrary truncation order would be 
more satisfactory. The recent formulation of Natarov et 
al. achieves just that. It consists in recasting the 2D- 
GLMT equation (53) in a block-type Fredholm form of 
the second-kind. One re-normalizes the beam shape co¬ 
efficients and the unknowns b^i using the formulas 

m 


^ni — Ji(A:o?’„) 

bnl — j 


(57) 


and solves the linear system for the new unknowns bni- 
Substituting in (52) yields 


C30 

bnl-SnlY, Y, 

n'^n I' —— 00 


Ji'{korn') I 
Mkorr,) 


SnlChnl- 


The modified matrix equation is therefore 

Tb = ao. 


(58) 


(59) 


with 


-*-nn' ~ '^ri 


(1 - 


Ji'(korn') 

Ji(korn) 


(60) 
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Figure 4. Transfer matrix for the triangle geometry shown in Fig. for kr\ = 5.3779. The truncation order of the matrix 
is chosen as Zmax = int[3fcri] + 1 = 17. (a) Magnitude of log|T| (b) Sub-matrix view of one off-diagonal block of log|T|. (c) 
Sub-matrix view of one off-diagonal block recast to a Fredholm form of the second-kind. 


The net effect of this re-normalization is that, for large 
1^1, the additional factor exponentially decays and com¬ 
pensates the exponential growth of the Hankel function 
with respect to that index. Conversely, for large \l\, the 
additional factor exponentially grows and compensates 
the exponential decay with respect to that index. Fol¬ 
lowing this rescaling, it can be shown that a truncated 
version of (59) converges to the exact solution of the 
scattering/eigenmode problem for ^max —t oo [6T]. In 
mathematical terms this stems from the fact that 


E . 

1,1'— — oo 


r^ll r 

nn' 


'nn' 


< oo, V (n, n). (61) 


In other words, the norm of the transfer matrix is finite 
regardless of the chosen truncation order ^max- Indeed, 
as can be seen from Fig. |^, the largest elements of off- 
diagonal blocks, transformed to a Fredholm form, are 
located near the diagonal of the blocks, whereas matrix 
elements remain bounded for large values of \l —I'\. This 
was not the case for the original system (53), as evidenced 
by Fig. |4 )d. 

For simplicity, we shall use the non-renormalized form 
of the central equation (53) of 2D-GLMT, unless indi¬ 


cated otherwise. Clearly, the theoretical discussion holds 
for the normalized version (59) as well. The computa¬ 


tions can always be carried out for T and b instead of T 
and b. Afterwards, one simply maps back the coefficients 
bni bni to reconstruct the electromagnetic fields. 


E. Scattering of arbitrary beams 


of linear equations (53). In fact, as stated by Gouesbet 
and Lock, the main difficulty behind Lorenz-Mie theo¬ 
ries is usually the computation of the beam-shape coef¬ 
ficients for a given excitation M- Although the calcu¬ 
lation of these coefficients is not our primary topic, it is 
an important ingredient for a complete implementation 
of the method. We therefore dedicate the next two sec¬ 
tions to the derivation of the beam-shape coefficients for 
two typical cases, namely a plane wave (section |lll E 1 ) 
and a focused beam (section III E 2) parametrized by the 
complex-source beam (CSB) technique. 


1. Beam-shape coefficients for a plane wave 

Consider the following incident field (plane wave) on 
the scatterer array 


LPiir) = (62) 

The incident wavevector k is defined by a modulus /cq = 
ky/eo and an angle of incidence 0. The position vector r 
is defined with respect to the origin of a global coordinate 
system. In the frame of reference of the n**' scatterer, one 
writes 


(^,(r) = = ^^ei^-Rr.^ikopAeos{e^-0)]^ (gg) 

and with the Jacobi-Anger expansion |62[ p. 361], one 
further has 

OO 

gzfeop„[cos(e„-e)] ^ ■ (64) 

1 — — 00 


The interaction of an incident beam - in principle of 
arbitrary shape - with an array of cylindrical scatterers 
can be readily modelled using 2D-GLMT. Given a set of 
beam-shape coefficients {a°;}, the scattered coefficients 
{bni} can be directly computed by solving the system 


The substitution of ( |64[ ) in (63) and comparison with 
(38) yields the beam-shape coefficients, identical to those 
found in [Tfl [551 [55] 


„0 _ ik R„ •/ —ii0 

- V’oe I e 


(65) 
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2. Beam-shape coefficients for a complex-source beam 


The basic form of a two-dimensional Gaussian beam 
(GB) propagating along the x axis, which satisfies the 
paraxial 2D Helmholtz equation, is given by 


Vg{x,y) = 


Trko{x - xr ) 


exp 




1 

2 X — ixB, 


( 66 ) 

The parameter xr is the Rayleigh distance of the GB and 
is proportional to the coherence length of the beam. The 
beam waist is located in the x = 0 plane. The physical 
appeal and usefulness of the Gaussian beam is well estab¬ 
lished. It is compatible with semi-analytical approaches 
such as ray-transfer (or ABCD) matrices and leads to a 
closed form solution of the Fresnel-Kirchhoff integral [65] . 
A Gaussian beam is also a good approximation of the ra¬ 
diation pattern of the fundamental mode of a rectangular 
waveguide [HI pp. 43-46] or of an optical fiber [67] • How¬ 
ever, the Gaussian beam solution suffers from two main 
drawbacks, a conceptual and a technical one. First, it 
is but an approximate solution of the Helmholtz equa¬ 
tion. One must use the paraxial approximation to obtain 
(66) from (12). Second, computing the beam-shape coef¬ 


ficients for the canonical Gaussian beam solution is not 
easily accomplished. 

To circumvent these two obstacles, one can use 
the complex-source beam (CSB) solution of the 2D 
Helmholtz equation as a closed-form incident field. This 
solution has been proposed in order to extend the valid¬ 
ity of the GB beyond the paraxial zone and exhibits the 


required cylindrical symmetry. Using the Green’s func¬ 
tion of the inhomogeneous Helmholtz equation for a point 
source located in complex space at coordinates x' = zxr 
and y' = 0, one obtains the CSB solution [65] 

(pi{r) = H^'''\kors), (67) 


with 

rs = [{y-y')^ + {x-x'fY/^ = [y^-b(x-'jxR)^]^/^. (68) 

This solution is continuous everywhere in the real plane 
except across the branch cut connecting the two singu¬ 
larities at (x,y) = (0,xr) and (x,y) = (0, —xr). The 
waist plane is thus located in the branch cut. It should 
be noted that to obtain a regular solution in that plane, 
one can combine linearly independent CSB solutions as 
described in [65] . 

The CSB solution reduces, up to a multiplicative con¬ 
stant, to a canonical Gaussian beam in the paraxial zone 
X ^ y. To show this, one can rewrite in the following 
way 


= ±(x - jxr)W1 -I- 


X ^ 0, 


(X - fXR)2 ’ 

and use the binomial approximation |69j to write 


± 


(X - fXR) -I- 


I 


y 


X ^ 0. 


(69) 


(70) 


2 X — ixR 

Moreover, using the asymptotic expansion of Hankel 
functions for large arguments, one can write 






7rfco(x - ixr) 


.. ( 1 2/^ 
exp iko X -I- - - 


2 X — ixr 


(71a) 


gT gfeoaJR 


X > 0, 


V>iix,y) 


7rfco(x - ixn) 


exp —iko ( X -I- - 


y 


2 X — fxR 


gTg koxn 


X < 0. 


(71b) 


Eq. (71) and Fig. [^clearly show the transverse Gaussian 
profile of the CSB in the paraxial zone. 

The next step in the derivation of the beam shape co¬ 
efficients for a focused beam is to expand the CSB on 
a basis of cylindrical waves centred on each individual 
scatterer. One rewrites (67) as 


<Pi{pn,0n) = H^'^\ko\r„ - 


(72) 


where r„ = {pn,0n) and rs„ is the vector pointing from 
the center of the n**' scatterer to the complex source 
point. We apply Graf’s addition theorem with the fol¬ 


lowing change of variables 

U = korsn, V = koPn, W = fcokn “ (73) 

where pn = |r„| and Xgn = |rsn| (Euclidean norms). In 
accordance with [60] . we write 

~ i^n T (74a) 

rsn = V(A„ - ixR)2 -I- Y^, (74b) 

thus 

^sn — i^n “f “t” ijjn “t” 4^)^y; 


(75a) 
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Figure 5. (Left) Phase fronts and Poynting vectors of the free-space CSB solution described by Eq. ( |67[ ). (Right) Amplitude 
contours of the CSB solution and comparison with the canonical Gaussian beam (dashed). The normalized Rayleigh distance 
is set to fcoXR = 9.65, the same value used in [2S]. Reproduced from [5S]. 


kn - rsnl = \/{Xn + “ IXr)'^ + {Vn + Y^), (75b) 


F. Eigenmode computations 


where (xmlJn) are Cartesian coordinates centred on the 
n*'' scatterer. Using ( [^ , one obtains 

cos(a„ + tt) = cos(6>„ -/r„), (76) 

where 

Y 

sinun = —, (77a) 

^ sn 

sin0„ = —. (77b) 

Pn 

Using these substitutions, the expansion takes the form 

OO 

^i{Pn,0n)= ^ i7,^^^(fcors„)Ji(fcoPn)e*'“", (78) 

1 — — 00 


cos Pn = 


Xn - iXK 


/I *^'1 

cos Un= - 

Pn 


and substituting an = On — Pn — tt, one finally obtains 
the beam-shape coefficients for a CSB 


a°, = (-l)'i7W(fcor,„)e-*'^" 


(79) 


where Vgn is given by ( 74b[ ) and /i„ by (77a). 

The expansion (78) is similar to Eq. (3) in reference 
with the difference that the condition 


< 


|U| is required in (SOI because the source is located inside 
a circular reflector, whereas (79) requires < |U|. 


In accordance with this condition, the convergence of (78) 


is limited to a disk not intersecting or touching the branch 
cut between (x,y) = (0,a;R) and (x,y) = (0, —xr). In 
other words, scatterers must not intersect or touch the 
branch cut for the expansion to hold. This condition is 
easily satisfied in any given setting. 


In the previous section, we have shown that the scat¬ 
tering wavefunctions of an array of cylinders can be com¬ 
puted by solving a linear system of equations. Eigenmode 
computations involve instead a homogeneous version of 
this linear system. Generally, this system is of the form 

T(A)b = 0, (80) 

i.e. the computation of lasing states is a non-linear eigen¬ 
value problem. One must compute the discrete set of 
(generally complex) eigenvalues A and eigenvectors b for 
which det(T) = 0. Since the non-linear eigenvalue is usu¬ 
ally a complex frequency or wavenumber, this problem 
amounts to finding the resonant frequencies associated 
to an infinite scattered amplitude in the presence of a fi¬ 
nite amplitude incident wave. Algorithms for solving the 
non-linear eigenvalue problem are described in IMZ2]. 

We now show that ( [^ can be readily adapted for the 
computation of the classical QB states of an array of 
cylinders, and of the CF states described in section [n B 2| 
The QB states of an array of dielectric scatterers satisfy 
the following Helmholtz equation 

-I- eok'^]p{pn, On) = 0, (Outside all cylinders) 

(81a) 

[V^ -I- enk'^]pipn, On) =0, < r„. (81b) 

Consequently, the matrix equation describing the quasi¬ 
bound modes is simply 

T(fc)b = 0, 


(82) 
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with the substitution kn •<— k^/e^ in (49). 

As stated in section II B[ QB states cannot accurately 
describe the steady-state lasing behavior of an array of 
active cylinders, even near threshold [S^]. This is due 
to the QB eigen-frequencies being complex everywhere 
outside the cylinders, resulting in exponential growth of 
the electromagnetic energy at infinity [^. To enable a 
more realistic treatment of eigenmodes, we have intro¬ 
duced in section IIB 2 a new kind of eigenstate central 
to SALT, the CF state. One of the appealing features 
of CF states is that they are readily computed using the 
Lorenz-Mie approach, under the single assumption that 
the cavity region is composed of a subset of the cylinder 
array. In other words, the medium surrounding all cylin¬ 
ders is passive, and some cylinders may also be passive. 
Accordingly, the wavevector is complex only inside active 
cylinders. The CF states therefore satisfy 


[V" + eok^]ip = 0, 
[V^ + enk‘^]ip = 0, 
[\/^ + e„Kikf]ip = 0, 


(Outside all cylinders) 

(83a) 

(Inside passive cylinders) 

(83b) 

(Inside the cavity region). 

(83c) 


The matrix equation describing the CF states is 


T(A:)b = 0, (84) 


with the substitution kn t— Ky/e^ in (49) if the n**' cylin¬ 
der is active, and t— k^/e^ otherwise. Note that, as 
always, the complex eigenfrequency K associated to the 
CF states depends on the value of the exterior frequency 

k. 


In summary, 2D-GLMT can be used to compute the 
eigenstates of arrays of coupled active cylinders, possibly 
all different in size and refractive indices. The method 
works equally well for the computation of the eigenfunc¬ 
tions and eigenfrequencies of QB states, meta-stable solu¬ 
tions of the 2D Helmholtz equation, and for the computa¬ 
tion of CF states. These “improved” eigenstates central 
to the SALT theory are more physically realistic solutions 
to the Helmholtz equation in the sense that they remain 
bounded at infinity, unlike QB states. Using CF states 
also enables the computation of resonances of photonic 
complexes in the case where some cylinders are active 
and other remain passive. This is incompatible with QB 
states computations, which amount to an active medium 
extending to the whole spatial domain. 


IV. NUMERICAL EXAMPLES 

This section is dedicated to numerical examples of 
both scattering and eigenmode computations using 2D- 
GLMT. It serves as a synthesis of all the theoretical con¬ 
cepts discussed in the previous two sections, including 
scattering of focused beams by coupled cylinders and 
eigenstates drawn from the SALT formulations. 


Two-dimensional photonic crystals (PhCs) based on 
arrangements of cylindrical scatterers are one of the 
prime candidates for modelling using 2D-CLMT [551 SZl 
1551 1751 Although PhCs are strictly speaking infi¬ 
nite periodic arrangements of scatterers, 2D-GLMT does 
not require a symmetric arrangement, nor is it meant 
to be applied to infinite structures. With these caveats 
in mind, we will sometimes refer to finite, periodic or 
non-periodic configurations as “PhC devices”. In section 
|IV A[ we present the numerical example of an integrated 
polarization filter as a representative scattering compu¬ 
tation using 2D-GLMT. In section |IV B[ the focus is on 
the more subtle eigenmode computations. These compu¬ 
tations are presented through the example of a so-called 
“photonic bandgap defect cavity”. 

Although our examples are restricted to devices com¬ 
prising a number of cylinders or cylindrical inclusions 
ranging from ~ 50 to several hundreds, they are by no 
means the only possible application of 2D-GLMT. Other 
notable instances of photonic and plasmonic devices that 
have been successfully modeled using 2D-GLMT include 
photonic molecules [Ml ESI ES] > integrated silicon lenses 
[261 El], nanowires [6l] and coupled-resonator optical 
waveguides m- 2D-GLMT can also be extended to ge¬ 
ometries characterized by an infinite extent in the z di¬ 
rection, such as micro-structured optical fibres (MOFs). 
The analytic formulation for MOFs can be found in [36] 
with details on the implementation and numerical exam¬ 
ples in m- In this “extended” formulation, the complex 
eigenfrequency fege of QB states described by equation 
(12) can be interpreted as a “transverse” propagation 
constant which can be used to compute a “longitudinal” 
propagation constant /3 using /3^ = fcg — fcgg, where fcg 
is the free-space wavevector. The modelling of MOFs 
constitute an important research topic, since these fibres 
can be used to accelerate charged particles using TM- 
polarized photonic bandgap defect modes analogous to 
those computed in section |IVB[ TM modes in hollow- 
core fibres are particularly appealing for electron accel¬ 
eration since they exhibit only a longitudinal electric field 
parallel to the fibre axis [38] . 


A. Scattering computations: polarization filters 


Our numerical application of a scattering computation 
is the characterization of integrated polarization filters 
based on engineered photonic lattices. This example 
serves two additional purposes. It illustrates the pos¬ 


sibility to use focused beams described by Eq. (67) in 


scattering computations. It also emphasizes the fact that 
the numerical approach described in section jTTTj allows to 
treat both orthogonal polarization components (TE and 
TM) independently. 

Our polarization filter is composed of cylindrical inclu¬ 
sions of refractive index n = I embedded in a medium of 
refractive index n = 2.76, corresponding to the effective 
value for a thin silicon slab. The radius of the inclusions 
















16 



0 2 4 6 8 10 0.0 0.5 1.0 

x/A 


Figure 6. Generation of a TM polarized Gaussian beam, (a) 
Optimized configuration (57 scatterers) and \E^\ field profile 
(arbitrary units). The target plane coincides with the upper 
limit of the x axis, (b) Gomparison of orthogonal polarization 
components along target plane. The energy of the polarized 
component of the beam is ~ 60 times higher than that of the 
unpolarized component. Configuration detailed in Ref. m¬ 


is set to r = 0.3A, where A is the lattice constant. This 
type of PhC geometry can be generically referred to as 
holes-in-slab (HIS). Although the filter is finite and non¬ 
periodic, its infinite counterpart - a square HIS lattice 
- exhibits a directional bandgap for the TE polarization 
around the normalized frequency k^K = 1.5. This im¬ 
plies that the TE component of the beam {Ey) is more 
strongly scattered, suggesting HIS based designs are suit¬ 
able for filtering this polarization out and favouring the 
TM polarization. This bandgap analysis, based on the 
plane wave expansion method, is detailed in m- 

Using optimization methods presented in [731 [74] , one 
is able to engineer a lattice configuration for simultaneous 
polarization filtering and beam shaping, that is obtaining 
a given beam shape after filtering. A configuration op¬ 
timized for the conversion of an incident Gaussian beam 
to a identical - but polarized - Gaussian beam is shown 
in Eig. [^ The configuration of cylindrical holes is shown 
on Eigj^, along with a false colour plot of the field am¬ 
plitude inside the device. In this example, the incident 
beam wavenumber is set to k^K = 1.76, chosen to fall in 
the vicinity of the directional bandgap of the infinite HIS 
lattice. Both the output and input beams’ half-width are 
set to Wo = 2.5A. The obtained output beam, shown in 
Fig. [§), deviates from a Gaussian amplitude profile by 
about 4% (RMS value), while the energy of the polar¬ 
ized component of the beam {Ez) is ^ 60 times higher 
than that of the non-polarized component {Ey). This 
corresponds to a degree of polarization of more than 98 
%. More details about the configurations, as well as the 
optimization of a TE polarizer, can be found in |74j . 

A few remarks concerning the scattering computation 
shown in Eig. [^ are in order. We did not use the renor¬ 


malized version of the matrix for scattering computations 
since it was not required for a numerically stable solution. 
Moreover, as stated in section [lHE2[ the complex-source 
beam approach precludes the positioning of scatterers in 
the beam waist. Therefore, we positioned all scatterers 
at x > r. For well-collimated beams, this is not a critical 
impediment. 

To illustrate the full potential of 2D-GLMT for pho¬ 
tonic structures characterization and optimization, it is 
instructive to consider the numerical cost of a computa¬ 
tion such as the one shown in Fig. [^ The dimension 
of the underlying transfer matrix is cr = A^(2Zi„ax + 1), 
where Zmax = 3int(fcr) -1-1 = 2. Since this polarization 
filter is composed of = 57 cylindrical scatterers, this 
translates to cr = 285. In other words, the computation 
of the scattered field shown in Fig. [^implies the solution 
of a square linear system of 285 equations, or equivalently 
the inversion of a 285 x 285 complex matrix. On a per¬ 
sonal computer, this operations takes a mere fraction of 
a second. Storing the larger transfer matrices associated 
to high Nkr values is also accessible, but may benefit 
from the use of parallel codes. Although this may seem 
restrictive, we stress the fact that large Nkr values do 
not lead to stability issues, especially considering the im¬ 
provement presented in section [IHD] which ensures that 
transfer matrices remain well-conditioned regardless of 
their size. 

The speed of the 2D-GLMT method is critical if the 
goal of a research project is the optimization of a pho¬ 
tonic device. Indeed, the fast characterization of cylin¬ 
der arrays enabled by 2D-GLMT allows it to be coupled 
with metaheuristic algorithms^ such as the genetic algo¬ 
rithm (GA). The combination of 2D-GLMT -|- GA has 
been successfully used in the past for the design of waveg¬ 
uide bends [S5], beam shapers [55] and integrated lenses 
[27] based on cylindrical scatterers. More generally, 2D- 
GLMT can be coupled with any optimization algorithm 
able to deal with a binary encoding of the problem. The 
scatterer configuration of Fig. was actually found us¬ 
ing a combination of 2D-GLMT and an algorithm called 
tabu search. More details can be found in two of our 
recent publications [71 [H]. Gontinuous optimization al¬ 
gorithms such as inverse harmony search m may also 
be combined with 2D-GLMT. For instance, one might 
decide to optimize not only the presence/absence of a 
scatterer, but also its geometric parameters such as its 
radius and refractive index. 


B. Eigenmode computations: photonic crystal 
cavities 

We now present numerical examples of GF states of 
coupled cylinder arrays Specifically, computing CF 


^ For all eigenmode computations presented in this tutorial, we 
use the renormalized version of the 2D-GLMT equations, as de- 
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Figure 7. (Top) Amplitude profile \ip\^ of a CF state of a 
photonic crystal defect cavity for an exterior frequency feA = 
1.885. The eigenvalue of this mode is AA = 1.885 — 0.0044i. 
(Bottom) Absolute value of the difference between the profile 
of the CF state and a neighbouring QB state with eigenfre- 
quency fcA = 1.885 — 0.0035i, or Q ~ 260 (log scale). 


photonic crystal laser. The bandgaps are forbidden fre¬ 
quency ranges arising from scattering and interference 
in a periodic structure. Introducing a defect in such a 
periodic structure may result in photons remaining con¬ 
fined inside, creating a so-called “photonic bandgap de¬ 
fect mode”, which is non-propagating. We focus our at¬ 
tention on point defects, as experimentally reported by 
Painter et al. [50] . 

The cavity geometry used is composed of 90 identical 
rods with a passive dielectric permittivity Cc = 13.18. 
The rods are arranged on a triangular lattice of which 
we consider only a hexagonal subset with a defect at its 
centre (see the overlay on Fig. [^. The radius of the rods 
is set to r = 0.3A, where A is the lattice spacing. The 
infinite counterpart of this finite photonic crystal array 
exhibits a complete bandgap for the TM polarization |0] . 

Before searching for the CF states of this cylinder ar¬ 
ray, it is useful to compute the QB states first since they 
can be associated one-to-one with CF states. It is subse¬ 
quently straightforward to look for CF states using the 
fact that K^{k') is usually located close to fege in the 
complex plane m- A false colour plot of a CF state 
found using this prescription is shown in Fig. One 
can see that the energy of the eigenmode is localized 
close to the point defect. The eigenvalue of this CF state 
{KA = 1.885—0.0044i) falls into the complete bandgap of 
the triangular rod lattice described earlier. On the bot¬ 
tom panel of Fig. we display the difference between 
the amplitude profile of this CF state and the associated 
QB state; it is smaller than 10“^ for the range consid¬ 
ered. Interestingly, one can see that the CF and QB 
states tend to differ mostly outside the contour contain¬ 
ing all scatterers. This is a testament to the fact that CF 
states remain bounded at infinity, unlike QB states, as 
also illustrated in Fig [l] In short, this numerical example 
shows that 2D-GLMT is versatile enough for the compu¬ 
tation of QB and CF eigenstates of arrays of cylinders. 


states using 2D-GLMT implies searching for eigenfre- 
quencies solutions of (84) and the associated eigen¬ 
functions using the scheme described in Appendix [Xj 
A countably infinite set of eigenvalue-eigenfunction pairs 
can be computed for every different value of the real ex¬ 
terior frequency A:, but only discrete values of k will yield 
a physically realistic TLM, as described in section |IIB| 
We shall not map CF states on individual TLMs, since in 
the case of uniform pumping this mapping is straightfor¬ 
ward. We refer interested readers to one of our previous 
articles for more details on this mapping m- 

A simple photonic crystal cavity is used as a repre¬ 
sentative cylinder geometry. Unlike a typical Fabry- 
Perot where confinement is achieved via total internal 
reflection, bandgaps are responsible for confinement in a 


scribed in section jmg This ensures a higher precision on the 
resonant frequencies and field patterns. 


C. Random lasers 

Random lasers 0 H] are non-conventional resonant 
structures where the usual laser cavity is replaced by a 
strongly scattering arrangement of particles embedded 
in an active medium, or of active particles in a host 
medium. These novel lasers sources are appealing struc¬ 
tures because of their potentially very small size, low- 
cost and ease of fabrication. The realization of random 
lasers using optically pumped active “powders” has been 
demonstrated experimentally mil]. It is also possible to 
fabricate disordered lasers in a more controlled way us¬ 
ing photonic-crystal-like structures. This was recently 
achieved via the etching of cylindrical inclusions in a 
GaAs planar waveguide, optically activated by layers of 
InAs quantum dots m- Several applications of random 
lasers have already been proposed, ranging from tumour 
detection [5] to optical trapping m- 

The recent interest in random lasers has given rise 
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to a number of modelling tools designed to harness the 
modal characteristics of strongly scattering random me¬ 
dia. SALT, for instance, provides an ideal framework 
for the study of complex lasing media, as a series of re¬ 
cent publications indicate [3M3ni[73[Hn]- Moreover, the 
potential of 2D-GLMT for the computation of modes of 
uniformly pumped random lasers was first proposed in 

[33]. 

We now consider the computation of the CF states of 
a large array of active cylinders via 2D-GLMT. We wish 
to illustrate the main numerical difficulties associated to 
the computation of GF and QB states of random lasers. 
Using Ref. [^ as a guideline, we assume that an array of 
identical, non-overlapping rods are uniformly distributed 
within a circle of radius Rout with a predefined surface 
filling factor. For the purpose of illustration, we consider 
one realization of an array with i?out = 40, normalized 
with respect to the radius of individual rods (r = 1). We 
set the filling factor to 0.2, for an array containing 320 
active scatterers with relative permittivity Cc embedded 
in air (e = 1). The geometry of the resulting array is 
shown in Fig. This arrangement of cylinders provides 
a useful, albeit simplified, model of a 2D random laser. 

The method for computing GF states described in Ap- 
pendix|^can be easily applied to the random array with¬ 
out significant modification. The main difficulty is rather 
numerical and lies in the much greater computation time 
for the determinant det[T] of the transfer matrix for a 
large number of cylinders. Also, the density of states 
in a given wavenumber interval is greater for a large ar¬ 
ray. Let us first set the permittivity of the cylinders to 
Ec = 4. Using a sufficiently fine sampling of the complex 
K plane (see Fig. [^), we are able to obtain estimates 
for the position of a number of GF states (see Fig. |^). 
These estimates can subsequently be refined using the 
procedure described in Appendix The field distribu¬ 
tion of one of these GF states is shown in Fig. For 
the parameters used (surface filling factor and refractive 
index), the scattering strength turns out to be relatively 
modest and the GF state intensity extends over a large 
fraction of the surface of the random laser. By increas¬ 
ing the scattering, for instance with a larger refractive 
index contrast between the cylinders and the surround¬ 
ing medium and/or a bigger filling factor, the eigenmodes 
of the system can become localized in an effective closed 
cavity. This spatial confinement is the principal attribute 
of an Anderson localized regime whereby the light inten¬ 
sity occupies a smaller, ultimately microscopic, region of 
the random array [311331 IMj- As an illustration, con¬ 
sider the result shown in Fig. [g which consists in a GF 
state of the geometry shown in Fig. [^ but for a value 
of Cc = 9. The Q factor of this resonance is approxi¬ 
mately twice that of the GF state shown in Fig. [^ and 
the intensity outside the cavity is smaller with respect to 
the peak intensity of the mode {Q = —Re[Ar]/2Im[Ar]). 
This is indicative of a transition towards a fully localized 
regime associated with minimal losses. For a thorough 
discussion of Anderson localization of light, we refer the 



Figure 8. A possible realization of a random laser composed 
of 320 identical active cylinders with a surface filling factor of 
0.2. Reproduced from |55| . 


interested reader to a recent topical review on the subject 

m- 

The computation of the field distribution with the res¬ 
olution shown in Fig. [^and[^is computationally inten¬ 
sive (several hours on a supercomputer) due to the large 
number of cylindrical harmonics that contribute to the 
total wavefunction at a given point in space (21niax + l for 
each scatterer). The transfer matrices associated to this 
geometry imply also important memory requirements. 
We should stress that eigenmode computations are al¬ 
ways more numerically costly than scattering computa¬ 
tions. The reason is clear: eigenmode computations often 
imply some kind of wavenumber sampling (as shown in 
Fig. 1^), giving rise to a much larger number of complex 
matrix operations ED- 

To conclude, a final point should be mentioned con¬ 
cerning a limitation of 2D-GLMT for the modelling of 
complex lasing media. Throughout our presentation, we 
have only considered photonic media composed of homo¬ 
geneous cylinders, i. e. piecewise constant refractive in¬ 
dex distributions. However, inhomogeneous and partially 
pumped resonators have been the object of many recent 
articles [54[ [79l [801 IS3H86] . The main conclusion of these 
publications, obtained using SALT, is that non-uniformly 
pumping a laser (random or not) is an ideal method to 
tune its emission properties. This partial pumping cor¬ 
responds to a dielectric distribution e(r) that varies in¬ 
side a single scatterer or resonator, possibly following a 
smooth function. Modelling inhomogeneous or partially 
pumped optical media represents an additional challenge 
since the spatial discretization of inhomogeneous regions 
is unavoidable. This unfortunately precludes the use of 
2D-GLMT as presented here. Despite this limitation, 
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Figure 9. Determination of CF states for the geometry shown 
in Fig. (a) Evolution of log | det[T]| in the complex K{k = 
2) plane, (b) Predictors of the position of various CF states 
of the random laser for £c = 4. These correspond to local 
minima in the surface defined by log | det[T]|. The solid dot 
indicates the position of the predictor of the CF state shown 
in Fig. 


2D-GLMT remains highly useful owing to its numeri¬ 
cal stability (especially considering the improvement pro¬ 
posed in section [III D ) and its relative speed in the ab¬ 
sence of spatial meshing. 


V. SUMMARY 

We have presented a detailed theoretical description of 
2D-GLMT, a numerical method dedicated to the mod¬ 
elling of scattering and eigenmodes of coupled cylinder 
structures. This method aims to exploit the internal sym¬ 
metry of multi-cylinder structures to allow the expansion 
of the electromagnetic field in a basis of cylindrical func- 



Figure 10. Field distribution of a single CF state with 
K = 1.9420 — 0.004i (Q ~ 240) for the geometry shown in 
Fig. E with €c = 4. The height coordinate is proportional to 
the intensity (TM polarization). 



Figure 11. Field distribution \ip\^ of a single CF state with 
K = 1.9184 — 0.002z {Q ~ 480) for the geometry shown in 
Fig.0 with Cc = 9. The height coordinate is proportional to 
the intensity (TM polarization). 


tions centred on every individual cylindrical scatterer. 
After a brief presentation of the basic differential equa¬ 
tions describing wave propagation in 2D passive (section 
0 and active (section |IIB| ) geometries, a detailed de¬ 
scription of the 2D-GLMT equations was presented (sec¬ 
tion the key part of which resides in the applica¬ 

tion of Graf’s addition theorem for cylindrical functions. 
We have also detailed how, after the application of appro¬ 
priate boundary conditions at the interface of individual 
cylinders (section IIIG), the 2D-GLMT equations can 


be reduced to matrix form. A recent improvement of the 
method consisting in reformulating the matrix equations 
into a Fredholm form of the second-kind form was also 
presented (section [ill D[ ). This improvement ensures that 
the transfer matrices can be truncated to an arbitrary 
order while still providing accurate numerical solutions. 

After this theoretical description, we have performed 
several numerical implementations to illustrate the po¬ 
tential of 2D-GLMT for the modelling of complex pho¬ 
tonic media, as well as the computational challenges as- 
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sociated to the use of the algorithm. These examples fall 
in two broad categories, namely scattering (section HIE) 
and eigenmode (section IIIF) computations. The char¬ 
acterization of PhC-inspired polarization filters was de¬ 
tailed as a repr esenta tive example of a scattering compu¬ 
tation (section IV A). This example also served to high¬ 
light the speed of 2D-GLMT, a useful feature when used 
in tandem with numerical optimization algorithms. This 
speed stems from the fact that 2D-GLMT does not in¬ 
volve a spatial discretization, or meshing, of the physi¬ 
cal space of the problem, unlike the more general FEM 
and FDTD algorithms. We have also presented a cele¬ 
brated application of 2D-GLMT, namely the computa¬ 
tion of modes of PhC cavities (section [lVB[ ). Finally, the 
computation of modes of random lasers was presented to 
push the algorithm to its numerical limits (section IV CI. 
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Appendix A: Computation of eigenmodes 


As stated in section [niF[ the computation of QB states 
is a non-linear eigenvalue problem (NLEP) of the form 

T(fc)b = 0, (AI) 


where T is a square complex matrix. A similar system 
for CF states is given by (84). NLEPs consists in find¬ 
ing the eigenvalues k and eigenvectors b satisfying ( Al ). 
The theory of linear systems prescribes that non-trivial 
solutions to this equation exist if and only if 


det[T(/fc)] = 0. (A2) 


This class of problems is termed “non-linear” since T de¬ 
pends non-linearly on k. While solvers are readily avail¬ 


able for linear eigenvalue problems, algorithms for solving 
the NLEP are more complicated. It is possible to solve 
the NLEP either using Newton’s method, or using lineari¬ 
sation methods that reduce the problem to a sequence of 
linear ones, as described in [701 [H]- 

In our work, we use Newton’s method to find the eigen¬ 
modes, as described in m The goal is to find roots of 
(A2) numerically. One iteration of Newton’s method can 
be written as 


det[T(fc„)] 

det[T(fc„)]'’ 


(A3) 


where k^ is an initial guess for the eigenvalue (complex 
wavenumber) and prime symbols indicate the derivative 
with respect to k. The derivative is computed using Ja¬ 
cobi’s identity [89] 

det[T]' = det[T]tr [T^^T'] , (A4) 

where T' is the element-wise derivative of T with respect 
to fc, and tr denotes the trace. This allows one to rewrite 
one step of Newton’s method as [72] 

" tr[T-i(fc„)T'(fc„)]- 


Newton’s method is known to converge quadratically to 
the roots of the determinant if one is sufficiently close to 
the said roots. To ensure convergence, the initial guesses 
for the eigenvalues fco must be chosen carefully. In our 
studies, we obtain these initial guesses by evaluating the 
value of det[T(fc)] in a region of interest of the complex 
k plane on a rectangular grid, as shown in Fig. for 
instance. Points of this grid corresponding to a change 
in sign of the function are subsequently used as initial 
guesses for Newton’s method, which in turn implies a 
numerical construction of T'. An alternative approach 
to computing det[T(A:)] is to use a singular value decom¬ 
position (SVD) [89|. Is consists in tracking the smallest 
singular value of T on a rectangular grid and using the 
minima of this “function” as initial guesses. This SVD 
approach can be useful when a geometrical symmetry of 
the cylinder array results in several degenerate non-linear 
eigenvalues, which can lead to a “false minimum” in the 
mapping of det[T(fc) ], as described in Ref. [57] . 

Once the roots of (A2) have been found using the New¬ 
ton algorithm, one simply uses a SVD routine to obtain 
the null eigenvector b, to reconstruct the electromag¬ 
netic field profile of the eigenmode using, for instance, 
Eq. (38b I. This method is general for QB and CF states 
and is not limited to 2D-GLMT computations. In fact, it 
can also be used in conjunction with other methods based 
on transfer matrices that imply NLEPs, for example the 
boundary element method (BEM) for cavity resonances 

[zniiTa. 
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